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Abstract 


The magnetometer has long been a reliable inexpensive sensor used in spacecraft momentum management 
and attitude estimation. Recent studies have shown an increased accuracy potential for magnetometer-only 
attitude estimation systems. Since the earth’s magnetic field is a function of time and position, and since 
time is known quite precisely, the differences between the computed and measured magnetic field 
components, as measured by the magnetometers throughout the entire spacecraft orbit, are a function of 
both the spacecraft trajectory and attitude errors. Therefore, these errors can be used to estimate both 
trajectory and attitude. Traditionally, satellite attitude and trajectory have been estimated with completely 
separate systems, using different measurement data. Recently, trajectory estimation for low earth orbit 
satellites was successfully demonstrated in ground software using only magnetometer data. This work 
proposes a single augmented Extended Kalman Filter (EKF) to simultaneously and autonomously estimate 
both spacecraft trajectory and attitude with data from a magnetometer and either dynamically determined 
rates or gyro- measured body rates. 


I. Introduction 

The magnetometer, due to its reliability and low cost, has been the focus of many studies in the recent past. 
Emphasis has been placed on using the magnetometer alone, without any additional input, to estimate the 
spacecraft trajectory (References 1, 2, and 3) and attitude (References 4 and 5). Studies have also been 
performed to determine the ultimate accuracy of the magnetometer in estimating attitude when accurate rate 
information is available (Reference 6). 

In using the magnetometer to estimate attitude, the spacecraft position is required to compute the reference 
magnetic field. In using the magnetometer to estimate position, including the spacecraft attitude improves 
the results. The data used to estimate either the position or the attitude is a function of the difference 
between the observed magnetometer measurements and the reference magnetic field. In this work we use 
this difference to estimate both the spacecraft attitude and position. This is an extension of the work 
performed by Shorshi and Bar-Itzhack (Reference 1) to add the attitude to the trajectory state vector. 

Many of the future missions, such as the Small and Mid-size Explorer Series and university class explorers, 
are looking for low cost and autonomous approaches to navigation and attitude estimation. This work could 
prove valuable to these missions as a prime navigation system, with coarse accuracy requirements, or a 
backup to a prime system where more stringent accuracy is required. 

In tiiis work we present the method of expanding the Extended Kalman Filter of Reference 1 to include the 
estimation of the spacecraft attitude, and the results of tests on the combined filter using simulated data. 
Incorporating the attitude into the filter requires an estimate of the rates. In this work we assume that the 
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rates would be provided by gyroscopes. A method similar to that of Challa (Reference 4) or Azor, Bar- 
Itzhack, and Harman (Reference 7) could be applied, though, in the absence of gyroscope data. 


II. Extended Kalman Filter Algorithm 


The EKF algorithm is based on the following assumed models: 


System Model: 

X =_f(X(t),t) + w(t) 

(1) 

Measurement Model: 

Zk= MXQtJ) +_v. 

(2) 


where w(t) is a zero mean white process, v* is a zero mean white sequence, and X(t) is the state vector 
defined as 


X T = [a, e, i, ft, co, 0, C d , g ] 

The first six elements of X(t) are the classical Keplarian elements which determine the spacecraft position 
and velocity, namely the semi-major axis (a), eccentricity (e), inclination (i), right ascension of the 
ascending node (ft), argument of perigee (co) , and true anomaly (0). C d is the drag coefficient and g 
represents the attitude quaternion. 

Measurement Update Stage: 

The linearization of equation (2) results in 


zjc = H k Xk + v k where H k = [Ho I HJ (3) 

Ho is the measurement matrix for the orbital states and is derived in Reference 1 , and H a is the measurement 
matrix for the attitude states. The derivation of H a is given in Appendix A. The effective measurement 
used by the filter is given as 


Zk ~ Em.k - B(X jctjJ 


(4) 


where B JRtk is the magnetic field vector measured by the magnetometer and B(Xk4k) is the estimated 
magnetic field vector as a function of the estimated state X k at time ^ . The dependence of B(X k ,t|J on the 
position and the attitude is seen in the derivation of equation (3) in Appendix A. 

The state update equation is 

X k (+) = XkW + KkZk (5) 

where is the Kalman gain computed according to 

Kk = P k (-)H k T [H k P k (-)H k T + RJ' 1 (6) 

R k is the measurement noise matrix and the covariance matrix is updated as usual with 

P k (+) = [I - K k H k ]P k (-)tI - KkHkf + KkRkKj (7) 

Equation (5) is used to update the orbital states, but not the attitude states. The update of the attitude states 
is done as follows. As shown above, the state vector contains the attitude represented by a quaternion. The 
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EKF estimates an error in the quaternion, expressed as a vector of three small angles, a, defined in 
Appendix A and derived in Reference 8. This error is combined with the current estimate of the quaternion 
to give an updated estimate of the quaternion, which is then propagated to the next time point. 

Propagation Stage: 

The filter dynamics model is given as 


where 


w 0 - 


A 0 

0 - 


where A k (X’ k ) = 



. W a_ 

0 

A a_ 


+i - A k (X’ k ) X’k+ 

X’ T = ja, e, i, £}, co, 0, C d , a ] 


( 8 ) 


Aq is the linearized transition matrix for the orbital states and is a function of the estimated orbital states, 
which are elements of X\ A* is defined in Reference 9. A a is the transition matrix for the attitude error, a, 
which are also included in X\ A a is based on the development from Reference 8. The transition matrices 
Ao and A a are first order approximations computed from the Jacobian F k (X v ) derived from the linearization 
of equation (1). 


The covariance matrix is propagated from time t k to time t k+1 using: 

Pw(-) = A k (XX+) k )P k (+)A k ( X (+) k ) T + Q k (10) 


Qic is the process noise covariance matrix for both the orbit and attitude states. The orbit states are 
propagated by solving equation (1) numerically without the noise component, as in Reference 8. The 
dynamics of the attitude states is linear. Assuming a constant angular velocity between gyro measurements, 
the attitude states are propagated using 

3ic.,<->=M k ( *> O') 


where 


and 
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and w = wi. 


( 12 ) 


T is the time between gyro measurements, w is the angular velocity vector, and the arguments 1, 2, 3 refer 
to the 3 components of w. Equations (10) through (12) are particularly suitable when testing with 
simulated data, because the rates are almost constant, with added noise. When the filter is applied to real 
data, equation (1) will be solved numerically without the noise component, as in Reference 8. 


III. Simulation 

A basic simulation was developed to test the EKF outlined above. The scenario consisted of simulating a 
spacecraft in low-earth-orbit with an earth-pointing attitude, i.e. maintaining a one revolution-per-orbit 
(RPO) attitude. The spacecraft axes, or body axes, are aligned with the orbital axes as defined in Figure 1. 
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Figure 1. Definition of Orbital Axes 


The rotation rate, w, used in equation (11), has a component only about the z<> axis. The derivation of the 
instantaneous rate experienced by the spacecraft is given in Appendix B. The axes labeled with T refer to 
the inertial coordinate system. Those marked with ‘o’ refer to the orbital coordinates. The quaternion 
represents the rotation from inertial coordinates to body coordinates. The attitude error, a, represents three 
small Euler angles around the body coordinates, which rotate the estimated quaternion to the true 
quaternion. The attitude displayed in the table below and in the results section is given in terms of Euler 
angles also. These Euler angles describe the attitude with respect to the orbital coordinates. Euler angles 
were chosen for the display since the true Euler angles are all zero. 

The parameters which define the baseline simulated orbit and attitude are given below. 


Parameter 

Truth 

A-priori 

estimate 

a (km) 

7000 

8000 

e 

0.05 

0.06 

i (deg) 

50 

54 

£2 (deg) 

90 

85 

co (deg) 

0 

5 

9 (deg) 

45 

50 

Q 

0.02 

1 

roll (deg) 

0 

10 

pitch (deg) 

0 

10 

yaw (deg) 

0 

10 
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IV. Results 


The simulation described above was run for 300,000 seconds. Noise was added to the simulated 
magnetometer data and to the simulated rate data. The magnetometer measurement noise was 2 milliGauss 
and the noise in the measured rate was 0.017 deg/sec . 

Figures 2 through 4 show the root-sum-square (RSS) error in the position estimate. The a-priori position 
error is 1453 km (computed from the orbital parameters given above). Figure 2 shows the error for the 
entire 300,000 seconds, approximately 5 1 revolutions (the orbital period is 97 minutes). Figure 3 shows the 
first 20,000 seconds. The error converges to less than 100 km within 10,000 seconds, which is roughly 1.7 
orbits. Figure 4 shows the final 50,000 seconds. The average converged position error is about 4 km. 




Figure 3. RSS Position Error - First 20,000 Seconds 
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Figure 4. RSS Position Error - Last 50,000 Seconds 

Figures 5 through 7 show the RSS attitude error. The attitude converges quickly, to less than 5 degrees, 
within 3,000 seconds as shown in Figures 5 and 6. Figure 7 shows that the average steady state error is less 
than 1 degree. 
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Figure 5. RSS Attitude Error 
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Figure 7. RSS Attitude Error - Last 50,000 Seconds 


Figures 8 through 10 show the RSS velocity error. The a-priori velocity error is 0.96 km/sec. Like the 
position, the velocity error converges within 10,000 seconds, as shown in Figure 8. Figure 9 shows that the 
error is less than 0.25 km/sec at the end of the first 5,000 seconds. Figure 10 shows that the steady state 
velocity error is approximately 0.004 km/sec. 
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Figure 8. RSS Velocity Error 



Figure 9. RSS Velocity Error - First 5000 Seconds 
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Figure 10. RSS Velocity Error - Last 50,000 Seconds 


Figure 1 1 shows the RSS measurement residuals for the first 70,000 seconds (the residuals are computed 
using equation (4) ). The average value is approximately 4 milliGauss. The residuals also converge quickly 
from an initial value of 1 86 milliGauss (RSS). 



Figure 11. RSS Residual 
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V. Conclusions and Future Work 


This work presents a single augmented Extended Kalman Filter that simultaneously estimates both 
spacecraft trajectory and attitude using data from magnetometers and gyroscopes. The results from the first 
test of this filter using simulated data, indicate that the filter can indeed estimate both the trajectory and 
attitude. Starting with errors (RSS) of over 1400 km in position and 10 degrees in attitude, the filter 
converged to less than 5 degrees in attitude within 3,000 seconds and to less than 100 km in position in 
10,000 seconds (1.7 orbits). The average steady state values are less than 1 degree for attitude and 4 km 
for position. The steady state velocity errors (RSS) are approximately 4 m/s ec and the average 
magnetometer residual is about 4 milliGauss (RSS). 

Further testing will be conducted both with simulated and real spacecraft data. The magnetic field varies 
more at higher inclinations. Therefore, the effect of the orbit inclination angle will be studied. Tests will be 
conducted as to the filter’s ability to estimate attitude and trajectory at low inclinations. The sensitivity to 
errors in Q will be examined. Shorshi and Bar-Itzhack (References 1 and 9) found that the estimation of £1 
was critical to the convergence of the position error. Additional errors will be introduced into the simulated 
data, e.g. magnetometer and gyro biases. The state vector will be expanded to include these biases and the 
ability of the filter to estimate these added states will be tested. Finally, tests with real spacecraft data from 
satellites such as the Gamma Ray Observatory, the Upper Atmospheric Research Satellite, and the Extreme 
Ultraviolet Explorer will be conducted. 
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APPENDIX A - Derivation of the Measurement Model 


The magnetic field vector can be resolved in the magnetic spherical coordinates, as shown in Figure A.l, as 
B t f = [B r , B eB , B«j> b ]* 



Figure A.l. Definition of the Magnetic Spherical Coordinates 

The magnetic field at the spacecraft location, computed using the IGRF magnetic field model and the 
estimated position, can be written as 

Y = D|,Df B F + n' (A.l) 

and the measured magnetic field vector, as measured by the magnetometer can be written as 

X m =D>DfB F + n m (A-2) 

where 

Dfc = the transformation from inertial to body coordinates 

p 

Dj = the transformation from magnetic spherical to inertial coordinates 
n' = the magnetic field model error 
n m = the magnetometer measurement error 

The effective measurement, z, is defined as follows 

2 = I m - Y = D^Bp + n m - D{,Df B F - X (A-3) 

Rewriting the transformation of Bp as 
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and 


DbDf B F = d£d, f B f + ACD^DfBp) 


(A.4) 


This leads to 


z = A(Dt,DfBp)+n (A.6) 

where 


A ( D b D f I F ) = AD^(Df B F ) + D^A(DfB F ) (A.7) 

The second term on the right hand side of equation (A.7) is the derivation of the measurement matrix for the 
orbital states given in Reference 9. The expansion of the first term leads to the measurement matrix for the 
attitude states. Rewriting that term as 


AD{,(DfB F ) = AD{,B I 


(A.8) 


where 


B, = the computed magnetic field vector in inertial coordinates 

The error in the transformation can be defined as the difference between the true body coordinates and an 
intermediate coordinate system, referred to as the computed body coordinate system. The matrix that is 

computed is £>{, , which is equivalent to a transformation to the computed body coordinate system, which 
can be written as 


D^-dJ-D^ (A.9) 

SO 

AD^D^-Dj, (A. 10) 

where D l b is the true transformation from inertial to body coordinates. For small attitude error we can 
assume that the matrix Dj? is composed of small angles, thus 


D 


b 

c 


therefore from equation (A. 11) 


0 - y 
y 0 

’•& <p 



-[a x] 


(A. 11) 


ADb = I - [a x] Dj, - D[=-(ax]D[ 

Substituting equation (A. 13) into the first term on the right-hand side of equation (A.7) yields 


(A. 12) 



AD b (D*B F ) = -[a x]D b B, =-[a x]B b =[B b x]a 


(A.13) 


Substituting equation (A.13) into equation (A.6) along with the measurement matrix for the orbital states, 
gives 


z= rB h x]a + HoX 0 + n_~ [Ho [Bb x]]x + n 


(A. 14) 


where Ho is the measurement matrix for the orbital states, Xo. and x is composed of both the orbital states 
and the small angular errors in the attitude, a. Since B b is not known, the magnetic field vector measured 
by the magnetometer is used instead. The combined measurement matrix is then given as 

H = [H 0 [B b x] ] = [H 0 H a ] (A. 15) 


APPENDIX B - Derivation- of Spacecraft Rotation Rate 


The instantaneous rotation rate about the spacecraft Zo is derived here from the orbital parameters which 



Figure B.l. Relationship Between Orbital Angles 

describe an elliptical orbit (the average rate is 1 RPO). Figure B.l defines the angles a, p, z, and J. The 
rotation rate, w z , is defined as 


The angle, a, can be written as 


but 


w 2 =bc 


a =0 + p 


(B.l) 


(B.2) 
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p = n / 2 + y 


(B.3) 


therefore 


oc = 0+ 7i/2+Y 


Equation (B.l) then becomes 


(B.4) 


w 2 =a -Q +y 


TTie relationship between the Keplarian elements, e, 6, and yis given as (Reference 10) 


then 


tan(y) 


e * sin(0) 

1 + e • cos(0) 


(B.5) 


(B.6) 


d d 


e ■ sin(0) 


: [1 + e ■ cos(0) J 

Performing the different.ation in equation (B.7) leads to the following equation 


w z =G ' 


cos 2 (y) 


1 + e cos(0) 


e 2 • sin 2 (0)O 
1 + e ■ cos(0) 


+ e • cos(0)0 


(B.7) 


(B.8) 


149 




